Self-Similar Spherical Collapse with Tidal Torque 
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N-body simulations have revealed a wealth of information about dark matter halos however their 
results are largely empirical. Using analytic means, we attempt to shed light on simulation results 
by generalizing the self-similar secondary infall model to include tidal torque. In this first of two 
papers, we describe our halo formation model and compare our results to empirical mass profiles 
inspired by N-body simulations. Each halo is determined by four parameters. One parameter sets 
the mass scale and the other three define how particles within a mass shell are torqued throughout 
evolution. We choose torque parameters motivated by tidal torque theory and N-body simulations 
and analytically calculate the structure of the halo in different radial regimes. We find that angular 
momentum plays an important role in determining the density profile at small radii. For cosmological 
initial conditions, the density profile on small scales is set by the time rate of change of the angular 
momentum of particles as well as the halo mass. On intermediate scales, however, p oc r -2 , while 
par" 



I. INTRODUCTION 

The structure of dark matter halos affects our under- 
standing of galaxy formation and evolution and has im- 
plications for dark matter detection. Progress in our un- 
derstanding of dark matter halos has been made both 
numerically and analytically. Analytic treatments began 
with work by Gunn and Gott; they analyzed how bound 
mass shells that accrete onto an initially collapsed object 
can explain the morphology of the Coma cluster 0, Q and 
elliptical galaxies Q. This continuous accretion process 
is known as secondary infall. 

Secondary infall introduces a characteristic length 
scale: the shell's turnaround radius r*. This is the ra- 
dius at which a particular mass shell first turns around. 
Since the average density is a decreasing function of dis- 
tance from the collapsed object, mass shells initially far- 
ther away will turnaround later. This characteristic scale 
should be expected since the radius of a mass shell, like 
the radius of a shock wave in the Sedov Taylor solution, 
can only depend on the initial energy of the shell, the 
background density, and time By imposing that the 
structure of the halo is self-similar - all quantities de- 
scribing the halo only depend on the background density, 
r ta (the current turnaround radius) , and lengths scaled to 
r ta - Bertschinger [ij and Fillmore and Goldreich (here- 
after referred to as FG) @ were able to relate the asymp- 
totic slope of the nonlinear density profile to the initial 
linear density perturbation. 

Assuming purely radial orbits, FG analytically showed 
that the slope v of the halo density distribution p oc r~ v 
falls in the range 2 < v < 2.25 for r/r ta *C 1. This devi- 
ates strongly from N-body simulations which find v < 1 
0,0 or v ~ 1-2 [1] at their innermost resolved radius and 
observations of Low Surface Brightness and spiral galax- 
ies which suggest y ~ .2 || and the presence of cores [Tol — 
\v\ . Though the treatment in FG assumes radial orbits 
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while orbits in simulations and observed galaxies contain 
tangential components, it is analytically tractable and 
does not suffer from resolution limits. Numerical dark 
matter simulations, on the other hand, do not make any 
simplifying assumptions and have finite dynamic range. 
Moreover, it is difficult to draw understanding from their 
analysis and computational resources limit the smallest 
resolvable radius, since smaller scales require more par- 
ticles and smaller time steps. It seems natural, then, to 
generalize the work done by FG in order to explain the 
features predicted in simulations and observed in galax- 
ies. This paper, in particular, investigates how non-radial 
motion affects the structure of dark matter halos. 

Numerous authors have investigated how angular mo- 
mentum affects the asymptotic density profile. Ry- 
den and Gunn analyzed the effects of non-radial motion 
caused by substructure [l3j while others have examined 
how an angular momentum, or a distribution of angular 
momenta, assigned to each mass shell at turnaround, af- 
fects the structure of the halo [T3 - |2ll | . Note that many of 
these authors do not impose self-similarity. Those that 
do assume that a shell's angular momentum remains con- 
stant after turnaround. 

This paper extends previous work by Nusser flij . 
Assuming self-similarity, he analytically calculated the 
structure of the halo in different radial regimes for shells 
with constant angular momentum after turnaround. He 
found that the inclusion of angular momentum allows 
< v < 2.25. According to Hoffman and Shaham [22j |. 
v depends on the effective primordial power spectral in- 
dex (d In P/d In k), which varies for different mass halos. 
For galactic size halos, Nusser's analytic work predicts 
v ~ 1.3, in disagreement with simulation results [f| @|. 
In order to address this discrepancy, we extend Nusser's 
work by including torque. We consistently keep track of a 
particle's angular momentum, allowing it to build up be- 
fore turnaround because of tidal interactions with neigh- 
boring protogalaxies [23] and to evolve after turnaround 
because of nonlinear effects within the halo. Moreover, 
we compare the predictions of our halo model to simula- 
tion results. 
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Self-similar secondary infall requires Q m = 1 since a 
nonvanishing fljy introduces an additional scale. Apply- 
ing self-similarity to halo formation in the ACDM model 
therefore requires approximations and a mapping to ha- 
los in an Einstein de-Sitter universe. We assume that 
the linear power spectrum and background matter den- 
sity p m today are equal in both universes so that the 
statistics, masses, and length scales of halos found in the 
two models are equivalent. Since the scale factor evolves 
differently in both universes, the halo assembly histories 
will differ. 

In section [TTl we define our self-similar system and 
torqueing parameters. In section IIII1 we set initial con- 
ditions and evolve the mass shells before turnaround. In 
section IIV[ we describe evolution after turnaround and 
then analyze the asymptotic behavior of the density pro- 
file at different scales in section [V] In section IVI1 we 
give numerical results, discuss the overall structure of the 
halo, and compare to N-body simulations. We conclude 
in section IVIII 



II. SELF-SIMILAR DEFINITIONS 

Here we explicitly define our self-similar system and 
derive constraints on the functional form of the mass dis- 
tribution within a halo and the angular momentum of 
particles in a particular shell. 

If the infall process is self-similar, then the halo's ap- 
pearance does not change once all lengths are scaled to 
the current turnaround radius. For our analysis, we de- 
fine the current turnaround radius as r ta (t) = Ct@ where 
both C and j3 are positive constants. The exponent /3, 
as we will find, depends on the initial perturbation spec- 
trum. 

The evolution of a particular mass shell must depend 
on time t and the shell's turnaround time t* . More explic- 
itly, assuming spherical symmetry, we have r = R(t, t*). 
We define a self-similar system as one in which every tra- 
jectory obeys the following scaling: 



R(At,AQ = A p R(t,U) 



(1) 



where A is a constant. The above implies that the tra- 
jectory of one mass shell with turnaround time t\ can 
be mapped to the trajectory of another mass shell with 
turnaround time t% = At\. The exponent /3 follows since 

J2(tl,tl)/fl(t 2 ,t2) = {tl/t 2 f. 

Each shell of a self-similar system must also follow the 
same equation of motion. From Newton's law, the radial 
equation of motion for a mass shell with angular momen- 
tum is given by: 



R(t, t* 



GM(R(t,Q,t) L 2 (R(t, f A t, t* 

' 1 ' ' (2) 



where dots denote derivatives with respect to the first 
argument, M is the mass of the halo interior to r and L 
is the angular momentum per unit mass of a particle in 
the shell. Note that we enforce the mass to not depend 
explicitly on t*, while the angular momentum can. As 
we will show below, this is physically motivated. From 
eq. (fT]), we find: 



R(At,AU) =A f} - 2 R(t,U) 



(3) 



Plugging in eqs. ([l) and (j3)) into eq. ([2]) and simplifying, 
we find: 



GM\A- p R(At,AQ,t 
R 2 (At,AU) 



R(At,At*) = - A 3/3 ~ 



L 2 ( A~PR(At, At*), t, U 

+ A 4 ^ 2 — i — ' (4) 

i? 3 (At,At*) 1 j 

Changing variables from R(t,U) to R(t,U)/CtP for the 
mass and angular momentum and rewriting eq. we 
find: 



GM(R(At,At*)/C(At)P,t 
R(At,At m ) = - A 3l3 ~ 2 v 



■ A 4/3_ 



R 2 (At, At*) 
L 2 (R(At 7 AU)/C{At)P,t,U 



R 3 (At,AU) 



(5) 



Relabeling coordinates and enforcing consistency with 
eq. (J2), we find the following constraints on the func- 
tional forms of the mass and angular momentum. 



M 



R 2 (t,U) 



i? 3 (t,t*) 



(R(t,U)/CtP,t) = A 3 P- 2 M(R(t,U)/Ct ,t/A) (6) 

L(R(t,U),t,t^ = A v - 1 L(iJ(t,t»)/Ct p s t/A ) t*/A) 

(7) 



With the above in mind, we define the angular mo- 
mentum per unit mass L of a particle in a shell at r, and 
the density p and mass M of the halo as follows. 



L{r,t)=B r Mlf{\t/U) (8) 

p(r,t) = PB (t)D(X) (9) 
4n 

M(r,t) = —p B (t)rl(t)M(X) (10) 

where A = r/r ta (t) is the radius scaled to the current 
turnaround radius and p B = l/6irGt 2 is the background 
density for an Einstein de-Sitter (flat Q m — 1) universe. 
Using eq. (fTJ) , it is straightforward to show: 



3 



III. BEFORE TURNAROUND 



\(t,AQ = \{t/A,U) (11) 

Eq. (TTT|) implies that if one can compute \(t,t*) for a 
particular mass shell at all times, then one also knows 
the position of all other mass shells, labeled by At* with 
varying A, at a particular time. This interpretation is 
very powerful and will be used later in order to calculate 
the mass profile after turnaround. 

If the mass profile M(X) also depended explicitly on 
£*, then the mass would not have to grow like the back- 
ground mass enclosed in the current turnaround radius. 
This is clearly not physical. Hence we suppressed the ex- 
plicit dependence on £*. On the other hand, we've kept 
the dependence on in the angular momentum in order 
to have this extra freedom. Inspired by tidal torque the- 
ory and numerical simulations, in eq. (|8]) we take / to 
be: 



The trajectory of the mass shell after turnaround de- 
termines the halo mass profile. In order to start inte- 
grating at turnaround, however, the enclosed mass of the 
halo must be known. For the case of purely radial orbits, 
the enclosed mass at turnaround can be analytically cal- 
culated [1, For the case of orbits that have a time 
varying angular momentum, we must numerically evolve 
both the trajectory and Ai(X) before turnaround in order 
to determine the enclosed mass at turnaround. 

The trajectory of a mass shell follows from Newton's 
law. We have: 

d 2 r _ GM(r,t) | L 2 (r,t) 
dt 2 r 2 r 3 

Rewriting eq. (|T3| in terms of A and £ = log(i/ii), where 
ti is the initial time, and plugging in eqs. ((8|), ([TO]) and 
(|T2"|). we find: 



/(A, t/t*) 



if t < U, 
if t > i*. 



(12) 



d 2 \ 



d\ 



The constant B sets the amplitude of the angular 
momentum at turnaround while 7 (tu) controls how 
quickly the angular momentum increases before (after) 
turnaround. Constraints on B, 7, and w will be dis- 
cussed in later sections. 

We've assumed that the halo is spherically symmetric. 
While simulated halos are triaxial [24j, the description 
above is meant to represent an average halo. Since there 
are no preferred directions in the universe, it should be 
expected that a statistically averaged halo is spherically 
symmetric. 

In the above, L represents the angular momentum per 
unit mass of all particles in the shell. We impose that 
all particles in the shell have orbital planes that are ran- 
domly distributed. This implies that the total vector 
angular momentum of the mass shell, and hence the to- 
tal angular momentum of the halo J, vanishes. Hence, 
while individual particles on a mass shell gain angular 
momentum in random directions throughout evolution, 
on average the mass shell remains spherical. Therefore, 
like we've assumed above, only one radial equation of 
motion is necessary to describe the evolution of the shell. 

Since our statistically averaged halo has a vanishing 
total angular momentum, this model cannot address the 
nonzero spin parameters observed in individually simu- 
lated halos [25|, [2(|. Nor can it reproduce the nonzero 
value o f (J 2 ) expected from cosmological perturbation 
theory [27H29|. However, j L 2 dm where dm is the mass 
of a shell, does not vanish for this model. We will use this 
quantity, which is a measure of the tangential dispersion 
in the halo, to constrain our torque parameters. 



— + (2/3-1)— +008- 1)A = -- 



2 M(X) 



9 A 2 



B 2 A -2 7 -3 

(14) 

The angular momentum before turnaround was chosen 
so that eq. (Tl4|) does not explicitly depend on £. This 
allows for a cleaner perturbative analysis. Since r is an 
approximate power law in t at early times, we still have 
the freedom to choose a particular torque model inspired 
by tidal torque theory. This will be discussed at the end 
of this section. 

In order to numerically solve the above equation, one 
must know A4(X), a function we do not have a priori. 
Before turnaround, however, the enclosed mass of a par- 
ticular shell remains constant throughout evolution since 
no shells cross. Taking advantage of this, we relate dr/dt 
to M. by taking a total derivative of eq. (ITUl) . We find: 



fi- - (3/3 - 2)Ct 



M 



M' 



(15) 



In the above, a prime represents a derivative taken with 
respect to A. Taking another derivative of the above with 
respect to time, plugging into eq. (|13|) and simplifying, 
we find an evolution equation for A4: 



(3/? - 2) 



/3(/3-l)A + (3/3-2)(/3 
M 2 M" 2M 



1 



M 



CM') 



9 A 2 



M 



(16) 



Given eqns. (114)) and (fl6|). we must now specify initial 
conditions when A > 1. We assume the following pertur- 
bative solutions for D(X) and A(£) valid at early times. 
M(X) follows from eq. (ITU|) . 
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FIG. 1. The variation of model parameter n with halo mass. 
Larger mass halos map to steeper initial density profiles. 



These represent the two solutions to the second order 
differential equation (eq. ITB"]) . For the first case, the 
turnaround radius grows faster than the Hubble flow 
while for the second it grows slower. Hence, the first 
solution represents the growing mode of the perturba- 
tion while the second is the decaying mode. Since we 
are interested in the growth of halos, we will only con- 
sider the growing mode from now on. Next, imposing 
p = 27 + 4 and enforcing equality between terms propor- 
tional to A 1_p in eq. (|16p. we find: 



9n 2 B 2 (p-3) 
2{p-n)(3n + 2p) 



(22) 



Comparing eqs. (fT5|) and (|2"2"|) . we see that the correc- 
tion to the initial mass caused by angular momentum is 
negative since p > n. This is expected since the angular 
momentum acts against gravity. 

Next, we find constraints for the parameters in eq. (|f 91) . 
Plugging in eq. (IT9l into eq. (fT4")l and setting terms linear 
in <5i, 82, Ai, and A2 equal to each other, we find: 



D{\) = l + 5 1 \- n + 5 2 \- p + ... (17) 

M(X)^X 3 (l + ^X- n + ^X-P + ..) (18) 
\ 3 — n 3 — p J 

A(0 = A e (2/3 -^(l + A ie Ql « + A 2 e Q2 « + ...) (19) 

In the above, n characterizes the first order correction 
to the background density. It is related to the FG pa- 
rameter e through n — 3e. It is also related to the ef- 
fective power spectral index n e ff = dlnP/dln/c through 
n = n c s + 3 [22J ■ Since n c s depends on scale and hence 
halo mass (Appendix [SJ , we have a relationship between 
n and halo mass. As Figure Q] shows, larger mass halos 
have larger n. This is expected since larger smoothing 
lengths imply steeper initial density profiles. As in FG, 
we restrict < n < 3 so that the density decreases with 
radius while the mass increases. We examine this whole 
range for completeness even though n > 1.4 corresponds 
to objects larger than galaxy clusters. The exponent p 
characterizes the next order correction to the background 
density caused by angular momentum. Consistency with 
our perturbative expansion (eqs.[T7]through[T5|) demands 
that we take n < p < 2n. However, it is straightforward 
to generalize to other cases (Section IIII A|) . As we show 
below, the constants {82, Ai, A2, «i, 012} are set by the 
equations of motion. The constants {Si, Ao} are set by 
boundary conditions. Plugging eq. (fT8|) into eq. (fl"6|) and 
enforcing equality between terms proportional to A 1- ™ 
we find two possible solutions. 

'HH) < 2o > 
>-§('-£) 



ai = - (23) 

Ai = V ( 24 ) 
11 — 6 

*2 = P \P - I J (25) 

A 2 = -^A/ (26) 
p o 

Eqns. (IT71) through (|2"6")l set the initial conditions for 
eqs (|T4"f and (lTo|) . We evolve A and M. and choose 
{Ao,<5i} such that turnaround occurs (dX/d^ = —A/3) 
when A = 1. The free parameters for evolution before 
turnaround are {n, B,p\. For the zero angular momen- 
tum case analyzed in FG, the enclosed mass is the same 
regardless of n; A-f (1) = (37r/4) 2 . Including torque, how- 
ever, we find that the enclosed mass depends on the pa- 
rameters B and p. Figure [5] shows a contour plot of 
lQAi(l)/9n 2 for n = 1. As expected, the enclosed mass 
at turnaround must be larger than the no-torque case 
in order to overcome the additional angular momentum 
barrier. B sets the amplitude of the angular momen- 
tum while p controls how the angular momentum grows 
in comparison to the mass perturbation. Larger B and 
smaller p correspond to stronger torques on the mass 
shell. Contour plots with different values of n give the 
same features. 

The above perturbative scheme ensures that the inclu- 
sion of angular momentum preserves cosmological initial 
conditions. Analyzing eq. (JT9"]) for £ —> —00, and using 
eqs. (|2"3")l and ([2"5jl . we see that since p > n, the angular 
momentum correction is subdominant to the density per- 
turbation correction. More importantly, at early times, 
the shell moves with the Hubble flow. Last, in order 
to be consistent with cosmological initial conditions, the 
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FIG. 2. Contour plot of 16M(1)/9tt 2 for n = 1 as a function of 
torquing parameters B and p. Smaller p and larger B result in 
larger torques on the mass shell. These larger torques require 
bigger enclosed masses at turnaround, in order to counteract 
the stronger angular momentum barrier. 



averaged halo, we can calculate (ct 2 ) using cosmological 
linear perturbation theory and compare to expectations 
from our model. 

Using the Zel'dovich approximation (30j . assuming a 
spherical Lagrangian volume with radius R, and working 
to first order, we find: 



(d*) M = Qa±&Mxl iax A\R) 



(28) 



where M is the mass of the halo, a is the scale factor, 
D is the linear growth factor, dots denote derivatives 
with respect to proper time, a; ma x is the lagrangian radius 
of the volume, R is the spherical top hat radius for a 
halo of mass M and A{R) is a time independent function 
defined in Appendix [B] which has units of length. Note 
that the scale factor a and D are the only quantities 
which vary with time. For a matter dominated universe, 
((a 2 )^ 1 / 2 oc t, just as in the White analysis. This is 
expected since the Lagrangian mass is time independent. 

Next we calculate eq. (|2T|) from the perspective of our 
model. Using eqns. ©, (|TU|) . and (fT2"j) and assuming first 
order corrections to M are negligible, we find: 



angular momentum of particles within a mass shell must 
vanish at early times. Imposing r oc i 2 / 3 , and plugging 
into eq. ©, we find that L cx i^+P/")/ 3 . Hence, for all 
values of p we consider, the angular momentum vanishes 
at early times and has a value of Brl/t* at turnaround. 
Note that shells which turn around later have larger an- 
gular momentum. 



A. Tidal Torque Theory 

According to Hoyle's tidal torque mechanism, mass 
shells before turnaround gain angular momentum 
through their interactions with the tidal fields of neigh- 
boring protogalaxies [23| . Peebles claimed that the angu- 
lar momentum of proto gala xies in an Einstein de-Sitter 
universe grows as i 5 / 3 [27] while Doroshkevich showed 
that for non-spherical regions, the angular momentum 
grows as t [28, 29]. White confirmed Doroshkevich's anal- 
ysis with N-body simulations [28j . Since the net angular 
momentum of our model's halo vanishes, we will instead 
use <7 2 , defined below, to constrain p and B. 



a 2 



dm\(r — t*o) x (v — Vq) 



(27) 



Vl 



In the above, we integrate over the Lagrangian volume Vl 
of the halo, r and v are the physical radius and velocity 
of particles within the halo, ro is the center of mass of the 
halo and Vq = v(tq). As described earlier in this section, 
A ^> 1 corresponds to early times when the halo is lin- 
ear. Therefore, since our model represents a statistically 



a 2 



B 



-in 



t 2 dr 
B 2 PB(t)rJ a (t) 



dr 



3 - 2 7 



t: 2 



27 ( 
min 



(29) 



The lower limit of integration sets an effective smoothing 
length which we choose to be A > 1 so that we only 
count shells that are still described by linear theory. The 
upper limit of integration is required since all the mass in 
the universe does not go into the halo. Since p = 2j + 4 
and n < p < 2n, then for the range of n we consider, 
27 < 3 and the angular momentum of the protogalaxy 
is dominated by shells close to r max . Equating eqs. (|28|) 
and (|29[) and assuming the first order corrections to r max 
in eq. (|19[) are negligible, we find p = 2n and: 



2n)M(l) {n - 1)/3 



A{R) 
R 



(30) 



Eqs. (HU) and ([3D) are derived in Appendix [Bj Note that 
the perturbative analysis, presented above, which is used 
to calculate A-^(l) is not valid for p = 2n. Redoing the 
analysis for this special case, we find: 



Q2 = - 



' 14 
9 

A, = —\ 



B 2 (2n-3) 



(7n- 17)(2n 



-2n 



B — - 



7(n~3) 2 



3) A2 



3 {n-Zf 



(31) 
(32) 

(33) 
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FIG. 3. The variation of model parameter B with halo mass. 
B sets the angular momentum of particles at turnaround. 



For the remainder of this paper we impose p — 2n, so 
that angular momentum grows in accordance with cos- 
mological perturbation theory, and set n according to 
the halo mass. Unfortunately, when comparing to N- 
body simulations (Section IVI Ap . eq. ([50)1 overestimates 
the angular momentum of particles at turnaround by a 
factor of 1.5 to 2.3. We discuss possible reasons for this 
discrepancy in Appendix IBl For convenience, -B1.5 (-B2.3) 
denotes B calculated using eq. (|3U| with the right hand 
side divided by 1.5 (2.3). As described above, M(l) in 
eq. (|30p depends on B. Therefore, in order to find B, we 
calculate B and M(l) iteratively until eq. (|30|) is satis- 
fied. 

The relationship between n and n e ff as well as tidal 
torque theory implies that {n,p, B} are all set by the 
halo mass. Figure [3] shows the variation of B with halo 
mass. A(R) increases with halo mass since more power at 
large scales is included; R also increases with halo mass. 
These two competing effects cause a slight variation in B 
over seven orders of magnitude in halo mass. 



IV. AFTER TURNAROUND 

Given the enclosed mass found at turnaround, we now 
solve for the trajectory and mass profile after turnaround. 
For convenience, we redefine the time variable to be 
£ = ln(i/i ta ), where t ta is the current turnaround time. 
The trajectory's evolution equation after turnaround, 
with the appropriate torque model feq. I12[) . is shown be- 
low. 



fx 
W 



Hl> 1)^- >-!)A 



2M(X) 
9 



B 2 



X 2 A 3 



; 2(w+l-20)£ 

(34) 



The torque model after turnaround was chosen to not 
explicitly depend on r, since r begins to oscillate on 
a much faster timescale than the growth of the halo. 
Nusser [l4[ and Sikivie et al. 17 [ focused on the case 



zu = 0. However, as was discussed before, this results in 
density profiles steeper than what is predicted by numer- 
ical simulations. 

There are a number of dynamical processes that can 
cause a particle's angular momentum to evolve after 
turnaround. Dynamical friction (3l| transfers the an- 
gular momentum of massive bound objects - like black 
holes, globular clusters and merging satellite galaxies - 
to the background halo. A massive black hole at the 
center of the halo that dominates the potential at small 
scales tends to make the velocity dispersion isotropic (32h - 

tBars [35l l3o| and supermassive black hole binaries 
[H| are also expected to perturb the dark matter 
velocity distribution. While the torque model proposed 
after turnaround is clearly very simplistic and may not 
accurately describe some of the above phenomena, it still 
allows us to get intuition for how torques acting on mass 
shells change the structure of the halo. 

Analytically calculating w is difficult since the halo af- 
ter turnaround is nonlinear. In Appendix [C] we show 
in a simplistic manner how w is sourced by substruc- 
ture and argue that dark matter dominated substructure 
should cause steeper density profiles than baryon domi- 
nated substructure. In order to properly constrain w, N- 
body simulations are required. This is beyond the scope 
of this work. 

The initial conditions for eq. (|34l) . enforced in the 
above section, are A(£ = 0) = 1 and dA/<i£(£ = 0) = — j3. 
As discussed before, self-similarity implies that all mass 
shells follow the same trajectory A(£). Hence, A(£) can 
either be interpreted as labeling the location of a partic- 
ular mass shell at different times, or labeling the location 
of all mass shells at a particular time. We take advan- 
tage of the second interpretation in order to calculate the 
mass profile. 

After turnaround, shells cross since dark matter is col- 
lisionless. Therefore, the mass interior to a particular 
shell does not stay constant. However, since A(£) spec- 
ifies the location of all mass shells at a particular time, 
the mass interior to a given scale is simply the sum of all 
mass shells interior to it. The mass profile is then given 
by|,l: 



2 r 00 

M{X) = -MCI) / d£exp[-(2/n)£]ff[A - A(£)l 
n Jo 

= M(l)^(-ir 1 exp[-(2/n)^] (35) 



where M(l) is the normalization constant found in the 
prior section, H[u] is the heaviside function, and £j is the 
ith root that satisfies A(£) = A. The above is straightfor- 
ward to interpret. The roots label shell crossings at a 
particular scale and the exponential factor accounts for 
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the mass difference between sheiis that turn around at 
different times. 

Since the trajectory and the mass profiie depend on 
each other, it is necessary to first assume a mass profile, 
then calculate the trajectory from eq. (f34|) and resulting 
mass profile from eq. Q35[) and repeat until convergence 
is reached. 

The density profile -D(A) is straightforward to derive 
using eqs. M EH and |35] We find 0,11 : 



D(X) = 



1 dM 

3A 2 dX 

2 M(l) 
3n A 2 



exp[-(2/n)&] 



d\\~ 



(36) 



ASYMPTOTIC BEHAVIOR 



Unlike N-body experiments, self-similar systems are 
not limited by resolution. One can analytically infer the 
asymptotic slope of the mass profile close to the origin. 
FG did this by taking advantage of adiabatic invariance, 
self-consistently calculating the mass profile, and ana- 
lyzing the limit of the mass profile as A — > 0. Below, 
we generalize their analysis to the case of particles with 
changing angular momentum. Unlike Nusser |Tj) , we do 
not restrict our analysis to the case zu = 0. 

Just as in FG, we start by parameterizing the halo 
mass and the variation of the apocenter distance r a . 



(37) 
(38) 



In the above r* is the turnaround radius of a mass shell 
which turns around at t*. It is possible to relate q and 
a to n by taking advantage of adiabatic invariance. The 
equation of motion for the mass shell is: 



M(r,t) = 


n(t)', 


r a _ 









dt 2 ~ GK ^ r + r 3 



(39) 



At late times, the orbital period is much smaller than 
the time scale for the mass and angular momentum to 
grow. Taking n(t) and L(t) to stay roughly constant over 
an orbit and integrating the above equation, we find the 
energy equation: 



drY 
dt) 



1 v a 

a — 1 



L 2 {t)(r- 



(40) 



The above relationship tells us how the pericenters r p 
evolve with time. Note that we only consider torquing 
models (eq. I12[) which give rise to bound orbits. This 



restriction on zu will be discussed below. Defining y 
fp/fa and evaluating the above at r 



r p , we find: 



l-y° 



i 



(«-i)£ 2 (t) 



(41) 



For zu > (<) 0, the angular momentum of particles in 
the mass shell increases (decreases). This gives rise to 
pericenters that increase (decrease). Hence at late times, 
the orbit of a mass shell with increasing angular momen- 
tum will circularize and have y ~ 1, while the orbit of 
a mass shell with decreasing angular momentum will be- 
come more radial, with y <C 1. With this in mind, we 
can now calculate the radial action in order to find how 
q relates to n and a. The radial action is given by: 



a — 1 



du[{l-u a - l )-A{y){u- 2 -l)] 1/2 (42) 

v(t) 

In the above we've assumed a > 1. Generalizing to 
the case a < 1 is straightforward. The special case a = 1 
will be addressed later. For y(t) <C 1, the above inte- 
gral is dominated by the region in which y(t) <C u <C 1. 
Over this region, the integrand is time independent and 
hence the same for all orbits. Therefore adiabatic in- 
variance implies n(t)r" +1 = const. For y(t) ~ 1, the 
orbit is circular, which implies the radial action vanishes 
and L 2 (t) = Gn{t)r% +1 . Using eq. (|3EJ), and noting that 
K(t) oc t s where s = 3/3 — 2 — a/3, we find at late times: 



_ j^{2^+^[a(l + n)-3}} if w > 
[a(l + n) - 3] if zu <c 



(43) 



3ra(a+l) 



For the specific case, zu < 0, taking advantage of y <C 1, 
the adiabatic invariance arguments above, and eqs. ([H]) 
and (|12p . we can rewrite eqn. (|41[) in the form y(t,t*) — 
yo(t/U) 1 , where: 



/ = 



vu if a > 1, 

2va/(a + 1) if a < 1. 



(44) 



and j/o r * is the pericenter of a mass shell at turnaround. 
Constant angular momentum after turnaround corre- 
sponds to zu = 0. This case was addressed analytically 
in fl4j | and numerically in |l7j |. 

We next take advantage of the functional form of the 
mass profile. Following FG, we define P(r/r ai y) to be 
the fraction of time a particle with apocenter distance r a 
and pericenter yr a , at a particular time t, spends inside 
r. 
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P(v,y) = (v<y) 

P( V ,y)= I ±iy\ (y<v<l) 

P(v,y) = l (v>l) (45) 

where 

{r v du if rv > 1 

Jy ((i- u »-i)-Afe)(«-2-i)) 1 /2 
r in if a < 1 

(46) 

We see that the presence of pericenters causes the new 
case v < y, which did not exist in the FG analysis. Self 
consistency demands that 

(_r_\ a = %t) [ M » dM^p ( r \ 

\nj M(r ta ,t) J M ta / 

(47) 

where M* is the mass internal to a shell that turns around 
at t* and Mta is the current turnaround mass. The in- 
tegral assigns a weight to each shell depending on how 
often that shell is below the scale r. Noting from eq. (fTU)) 
that 



tt a only contribute. As a result, y(t,t*) <C 1. Using eq. 
(gU), we then find: 

y(t,Q = y (Pl =yo( J —) (si) 

where 8 = l/(q — /3). For bound mass shells, q — (3 < 0. 
Therefore, since 8 > 0, the first argument of P in eq. (|49)l 
increases while the second decreases as we sum over shells 
that have turned around at earlier and earlier times (u — > 
oo). For r/r ta "C yo, mass shells which most recently 
turned around do not contribute to the mass inside r / r ta 
since we are probing scales below their pericenters. Mass 
shells only begin to contribute when the two argument of 
P are roughly equal to each other. This occurs around: 

u = 2/1 = [yo {r/r ta ) ) (52) 

Hence, we can replace the lower limit of integration in 
eq. (14T)|) with j/i . We next want to calculate the behavior 
of eq. (|49p close to y\ in order to determine whether the 
integrand is dominated by mass shells around y\ or mass 
shells that have turned around at much earlier times. 
The first step is to calculate the behavior of P(u, y) for 
u w y. We find: 



M. = M ta 



3/3-2 



1/2 



(48) P(u, y) oc u 1/2 (l - y/u) 1 ' 2 x \ ^_ 



if a > 1, 
"/ 2 if a < 1. 



(53) 



using eq. (|38|) and transforming integration variables, we 
find: 



where 



a— k 
— ) = A' 



k = 



r/r tB . 



du 

,,1 + fe 



P(u,y(t,Q) (49) 



2 + n(2-3<?) 



(50) 



As u increases, the above integral sums over shells with 
smaller f*. Since the pericenter of a shell evolves with 
time, the second argument of P depends on u. The de- 
pendence, as we showed, varies with torque model (sign 
of vj); hence we've kept the dependence on u implicit. 
Next we analyze the above for certain regimes of r/r ta , 
and certain torquing models, in order to constrain the 
relationship between a and k. 



A. r/rta < yo, H7 < 

For vj < 0, particles lose angular momentum over time. 
When probing scales r/r ta -C yo, mass shells with t* <C 



Given the above, we evaluate the indefinite integral in 
eq. (l49|) . noting that y is a function of u (eq. [5Tj) . For 
m ~ 7/1, we find: 



P 




.i-fc 



if a > 1, 



oc (u/ yi - If' 2 \%_ k _ a/2 

1 2/x if a < 1. 

(54) 

Now comes the heart of the argument. Following the 
logic in FG, if we keep ujy\ fixed and the integrand blows 
up as r/r ta — > 0, then the left hand side of eq. (|4*9"f must 
diverge in the same way as the right hand side shown in 
eq. (|54l) . Therefore, using eq. ([52"]) : 



-k = 



5(1 -*)/(! + *) 



if a > 1, 



5(3/2- fc -a/2) (1 + 5) ifa<l. 



(55) 



Otherwise, if the right hand side converges, then the inte- 
grand will not depend on r/r ta as r/r ta — ► 0. Therefore, 
the left hand side cannot depend on r/r ta either, which 
implies a — k. Solving the above system of equations 
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for a given eqs. (l43l and (1501 and making sure the solu- 
tion is consistent, (ie: using eq. (|55[) only if the integrand 
diverges), we find: 

For n < 2 : 

1 + n - -v/Cl + ") 2 + 9nnj(nn7 - 2) 

a = 

3nn7 



fc = 



<7 = 



1 + n + 3mu — -\/(l + ?i) 2 + 9nw(nzu — 2) 
ntc7(4 + n) 

1 + n — 3ntu — v0~+ n ) 2 + 9nn7(ntn — 2) 



3n 



For n > 2 : 



1 + n 



, <z = 



(56) 



The above solutions are continuous at n = 2. Moreover, 
taking the no-torque limit (w — > 0) for n < 2 gives the 
same solutions as n > 2, which is consistent with analytic 
and numeric results from Taking the limit, w — > 

—oo reproduces the FG solution, as expected, since the 
shell loses its angular momentum instantly. The solution 
for n > 2 is independent of 137. This is because the mass is 
dominated by shells with turnaround time t* <C U a which 
have effectively no angular momentum. In other words, 
for w < 0, the solution should only depend on torquing 
parameters when the mass is dominated by shells that 
have turned around recently. 



B. 



r/r ta < yo, no > 



1 + n — 3nro 



For w > 0, the angular momentum of particles increase 
with time. As mentioned above, when probing scales 
f/fta *C 2/0: mass shells with t» <C t ta only contribute. 
As a result, y(t,t*) ^ 1. In other words, the orbits are 
roughly circular. We can therefore replace the lower limit 
of integration in eq. (|49[) with 1 since mass shells will only 
start contributing to the sum when u ~ y ~ 1. Hence, 
the right hand side of eq. (j49|) does not depend on r/rta,, 
which implies a — k. Using eq. (1431) and (j50")l , we find: 



r = 2m , for < n < 3 
(57) 

The no torque case, 137 = 0, is consistent with the 
analysis in the prior subsection. The singularity zu — 
(1 + n) /3n implies q = ft. This physically corresponds to 
the orbital radius of mass shells increasing at the same 
rate as the turnaround radius and results in orbits that 
are not bound and a cored profile where there are no 
particles internal to a particular radius. This breaks the 
assumption of a power law mass profile (eq. I37p ; hence 
we only consider w < (1 + n)/3n. 

Unlike the Nusser solution, certain parameters give 
a > 3, which corresponds to a density profile which con- 
verges as r — > 0. Since the angular momentum acts like a 
heat source dp/ dr > is dynamically stable and physical. 



C. y < r/n a < 1 

In this regime, we are probing scales larger than the 
pericenters of the most recently turned around mass 
shells. As a result, P{u, y) is dominated by the contribu- 
tion from the integrand when n» y. Therefore: 



P(uv) ocT if«>l. 
Hence the integral in eq. (|4^]) becomes: 



(58) 



du 1 u 



i-fc 



,i+fc 



if a > 1, 

3/2-k-a/i ifa< l. 



(59) 



Following the logic in the prior section, if the integral 
diverges as r/r ta — > 0, then we set the exponents on the 
left hand side and right hand side equal to each other 
so that both sides diverge in the same way. If the inte- 
gral converges, then the left hand side cannot depend on 
r/rtu, which implies a — k. Given these arguments, and 
imposing consistency with the above inequalities on a to 
find the appropriate ranges for n, we find: 



6 n - 2 
a = 1 , k = — — , q = — , for n < 2 



4 + n 
3 

a = k = — — , (? = 0, 
1 + n 



3ri 



for n > 2 



(60) 



The above is exactly the FG solution. We expect to 
recover these solutions since we are probing scales larger 
than the pericenters of the most massive shells, where 
the angular momentum does not affect the dynamics. 

This section assumed a / 1 and yet, for certain parts 
of parameters space, eqs. (|56|) . (|57| and (|60| give a = 1. 
However, since the solutions are continuous as a — >• 1 
from the left and right, then the results hold for a = 1 
as well. 



VI. STRUCTURE OF THE HALO 

In this section, we discuss the radial structure of galac- 
tic size halos and compare directly to numerical N-body 
simulations. Note however, that the mass of a halo is not 
well defined when our model is applied to cosmological 
structure formation since it is unclear how the spherical 
top hat mass which characterizes the halo when it is lin- 
ear relates to the virial mass which characterizes the halo 
when it is nonlinear. For halos today with galactic size 
virial masses, we assume the model parameter n which 
characterizes the initial density field, is set by a spherical 
top hat mass of 10 12 M Q . As described in prior sections, 
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specifying the top hat mass also sets model parameters B 
and p. Before comparing directly to N-body simulations, 
we first describe how zu influences the halo. 

Figure 0] shows the mass M{X) and density profiles 
D(X) for galactic size halos {n — 0.77, p — 2n, B1.5 = 
0.39} with varying vj. The spikes in the density pro- 
file are caustics which form at the shell's turning points. 
They form because of unphysical initial conditions; we as- 
sume each shell has zero radial velocity dispersion. The 
structure of the halos naturally break down into three 
different regions. The dividing points between these re- 
gions are roughly the virial radius (r v ) and yo^ta, the 
pericenter of the mass shell which most recently turned 
around. 

As in dark matter N-body simulations, we associate 
the virial radius with r2oo? the radius at which M(r2oo) = 
8007rpsr| 00 /3 is satisfied. Numerically, we find that the 
virial radius occurs near the first caustic (A ~ .18). For 
r > r v , the mass profile flattens and then starts to in- 
crease. The flattening is equivalent to what is seen on 
large scales in N-body simulations where p oc r~ 3 . The 
mass profile then starts to increase again because at large 
radii, where A> 1, the density is roughly constant, which 
implies M. oc A 3 . For r ~ r v , it is difficult to make an- 
alytic predictions for the mass profile because adiabatic 
invariance breaks down. In other words, the mass of the 
halo and angular momentum of a shell change on the 
same time scale as the shell's orbital period. 

As discussed in the prior section, for yor ta <C r <C r v , 
we can take advantage of adiabatic invariance to infer the 
logarithmic slope of the mass profile. Since this regime 
probes a scale much larger than the pericenters of the 
mass shells, the angular momentum does not affect the 
dynamics and we recover the FG solution. For our par- 
ticular choice of n = 0.77, this gives an isothermal profile 
with p oc r~ 2 . However, since n < 2 for all collapsed 
objects today (Figure [1]) and based on the results of FG, 
our model predict that all halos are isothermal in this 
regime. 

Last, for r/r ta <c yo, angular momentum begins to 
play a role and the halo starts to exhibit different features 
than the FG solution. The behavior is very intuitive. The 
mass of a particular shell does not contribute to the in- 
ternal mass when probing radii less than the pericenter 
of that mass shell. Therefore, as one probes radii smaller 
than the pericenter of the most recently turned around 
mass shell, one expects a steeper fall off than the FG so- 
lution, since less mass is enclosed interior to that radius. 
Moreover, varying w varies the pericenter of mass shells 
over time. Increasing (decreasing) angular momentum, 
vj > (<) 0, causes the pericenters to increase (decrease) 
over time. This results in profiles which are steeper (shal- 
lower) than the no-torque case. 

It is informative to find how the transition radius yo 
depends on model parameters. Since the mass and angu- 
lar momentum grow significantly before the first pericen- 
ter, we can only approximately determine this relation- 
ship. We assume that the profile is isothermal on large 
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FIG. 4. The mass and density profiles for galactic size halos 
{n = 0.77, p — 2n,B\.5 — 0.39} with varying zu. The value 
of zu changes how pericenters evolve with time and thereby 
affects how many shells at a particular scale contribute to 
the internal mass. The above numerically computed profiles 
match analytic predictions. The virial radius (r v ) and first 
pericenter passage (yorta) are labeled for clarity. 



scales and the halo mass and shell angular momentum 
are fixed to their turnaround values. For yo <SC 1 , the 
transitional radius yo solves the transcendental equation 
ygln(yo) = -9B 2 /4M(1). As expected, B -> repro- 
duces yo — >■ 0. As shown in Figure [21 -M(l) varies with 
{B,p}. However, for reasonable parameter values, the 
mass normalization is changed at most by a factor of 2. 
Therefore, yo most strongly depends on B, and p has a 
negligible effect on the structure of the halo. As seen in 
figure 21 yo should also depends on vj. The above ap- 
proximation neglects this dependence since we assumed 
the angular momentum is set to the turnaround value. 
Figure \5\ shows the radius of a mass shell, normal- 
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FIG. 5. The radius of a mass shell, normalized to its 
turnaround radius, for a galactic size halo {n — 0.77, p = 
2n, Bi.5 = 0.39}, plotted vs time. The top panel shows a 
shell with particles that have decreasing angular momentum 
(zu = —0.2). The middle panel shown a shell with particles 
that have constant angular momentum (w = 0). The bottom 
panel shows a shell with particles that have increasing angular 
momentum (zu — 0.2). 
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FIG. 6. A phase space diagram of a galactic size halo {n = 
0.77, p = 2n, Bi.5 = 0.39} at the current turnaround time. 
Velocities are normalized to the turnaround time and radius 
of each shell. The top panel shows a halo with zu = —0.2. The 
middle panel shows a halo with zu = 0. The bottom panel 
shows a halo with vj = 0.2. 
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ized to its turnaround radius r*, for a galactic size halo 
{n = 0.77, p = 2n,B 15 — 0.39}, as a function of time. 
In the top panel, particles in the shell lose angular mo- 
mentum (zu = —0.2), in the middle panel the angular 
momentum remains constant (zu — 0), while in the bot- 
tom panel particles in the shell gains angular momentum 
(zu = 0.2). As expected, the pericenters in the top panel 
decrease with time while the pericenters grow in the bot- 
tom panel. Hence, the orbits of particles with decreasing 
angular momentum become more radial while those with 
increasing angular momentum become more circular. 

Notice that the period of oscillation also varies for dif- 
ferent zu. The period of oscillation is set by the shell's 
apocenter r a and the mass internal to r a . Using the adi- 
abatic invariance relations we found in Section [V] and 
assuming Kepler's third law, we find that the period of 
the orbit P oc r\ for shells with decreasing angular mo- 
mentum and P oc r a for shells with increasing angular 
momentum. Moreover, from eqs. (|56[) and (1601) . r a de- 
creases (increases) with time for zu < (>) 0. Though 
Kepler's third law doesn't hold for this system, it still 
gives intuition for the above results. 

Figure [Fj] shows the phase space diagram for a galactic 
size halo {n = 0.77, p = 2n,B 1 . 5 = 0.39}. In the top 
panel, the particles in the shell lose angular momentum 
(zd = —0.2), in the middle panel the angular momentum 
remains constant (zu = 0), while in the bottom panel the 
particles in the shell gain angular momentum (zu = 0.2). 
The diagram labels the phase space point of every shell at 
the current turnaround time. All radial velocities are nor- 
malized to the shell's turnaround time and radius r*. 
Unlike FG, the presence of angular momentum results in 
caustics associated with pericenters as well, which can be 
seen in the lower panel of Figure @] In addition, since an 
increasing angular momentum results in increasing peri- 
centers, the pericenter caustics are more closely spaced in 
the lower panel than in the upper panel. Moreover, the 
amplitude of the radial velocity is smaller in the lower 
panel because orbits are circularizing. The phase space 
curve appears to intersect itself because we did not plot 
the tangential velocity component. In full generality, the 
distribution in the phase space (r,v r ,Vt = L/r) for our 
model is a non-self-intersecting one-dimensional curve. 



A. Comparing with N-body Simulations 

In this subsection, we compare the density profile of 
our model's halo to empirical fits inspired by N-body 
simulations. We first numerically calculate the density 
profile for a galactic size halo with zu = 0.12. This 
value of zu was chosen so that p oc r _1 on small scales. 
We then compute the spherically averaged density in 50 
spherical shells equally spaced in log 10 r over the range 
1.5 x 10~ 4 < r/r v < 3, and take r v = r 2 oo (defined 
above). This is the same procedure followed with the 
recent Aquarius simulation Q. Next, we calculate r_2, 
the radius where r 2 p reaches a maximum. For our halo, 
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FIG. 7. Spherically averaged density profile for the secondary 
infall model compared with NFW and Einasto profiles. The 
secondary infall model is calculated for a galactic size halo 
{n = 0.77, p — 2n} with vj — .1167. In the top panel we chose 
Bi.5 — 0.39 while in the bottom panel we chose B2.3 = 0.26. 



as discussed above, the profile is isothermal over a range 
of r. Moreover, the maximum peaks associated with the 
caustics are unphysical. So, we choose a value of r_ 2 
in the isothermal regime that gives good agreement with 
the empirical fits. Changing r_2 does not change our 
interpretation of the results. 

In Figure [7] we compare our spherically averaged den- 
sity profile to NFW and Einasto profiles. We plot r 2 p in 
order to highlight differences. The NFW profile is given 



p(r) 



4p_ 



(r/r- 2 )(l + r/r- 2 ) 2 
while the Einasto profile is given by: 



(61) 
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-1.42 



In p(r)/p_ 2 = (-2/a B )[(r/r_ 2 ) < « 



(62) 



where p-i is the density of our halo at r_2 and ag, known 
as the shape parameter, sets the width of the r 2 p peak. 
In the top panel, we use B\.§ — .39 while in the bottom 
panel, we use -B2.3 = -26. We choose cxe = 0.159 since 
this value was used in Figure 3 of Navarro et. al. [6j. 

We see that the secondary infall model works surpris- 
ingly well. The peaks are a result of the caustics that 
arise because of cold radial initial conditions. The first 
spike on the right comes from the first apocenter passage 
while the second comes from the first pericenter passage. 
The location of pericenter is most strongly influenced by 
the model parameter B. Hence, the isothermal region is 
smaller in the top panel than in the bottom panel since 
particles have less angular momentum at turnaround in 
the lower panel than in the upper panel. The parameter 
B then plays the same role as a^; it sets the width of 
the isothermal region. If we assume N-body simulations 
faithfully represent dark matter halos, then Figure [7] im- 
plies that our estimate of B in eq. (|3"0")) overestimates the 
actual value by 1.5 to 2.3. We discuss possible reasons 
for this in Appendix [B] 



VII. DISCUSSION 

N-body simulations reveal a wealth of information 
about dark matter halos. Older simulations predict den- 
sity profiles that are well approximated by an NFW pro- 
file [39(, while more recent simulations find density pro- 
files that fit better with a modified NFW profile [40j or 
the Einasto profile Q. In an attempt to gain intuition 
for these empirical profiles, we've generalized the self- 
similar secondary infall model to include torque. This 
model doesn't suffer from resolution limits and is much 
less computationally expensive than a full N-body simu- 
lation. Moreover, it is analytically tractable. Using this 
model, we were able to analytically calculate the density 
profile for r/r ta <C yo an d yo <C r/r ta <C 1. Note that 
the self-similar framework we've extended predicts power 
law mass profiles on small scales. Hence, it is inconsistent 
with an Einasto profile. 

It is clear from our analysis that angular momentum 
plays an essential role in determining the structure of 
the halo in two important ways. First, the amount of 
angular momentum at turnaround (B) sets the width of 
the isothermal region. Second, the presence of pericenters 
softens the inner density slope relative to the FG solution 
because less mass shells contribute to the enclosed mass. 
Moreover, the interior density profile is sensitive to the 
way in which particles are torqued after turnaround (w) . 

If we assume that w is constant for all halos, then this 
secondary infall model predicts steeper interior density 
profiles for larger mass halos. More specifically, if we use 
the value of w — 0.12 which gave p oc r^ 1 for galactic size 



halos, then p oc r -t bb for a 10 Mq halo and p oc r~ 
for a 10 15 Af Q halo. This trend towards steeper interior 
slopes for larger mass halos, and hence non- universality, 
has been noticed in recent numerical simulation s I4ll . |42( 
as well as more general secondary infall models [43j. On 
the other hand, if we assume that all halos have p oc 
as r — > 0, then w must vary with halo mass. More specif- 
ically, halos with mass M < 10 9 M Q must have particles 
which lose angular momentum over time (w < 0) while 
halos with mass M > 1O 9 M must have particles which 
gain angular momentum over time (w > 0). In other 
words, in order for our self-similar framework to predict 
universal density profiles, w must conspire to erase any 
dependence on initial conditions. A more thorough treat- 
ment requires the use of N-body simulations, which is 
beyond the scope of this paper. 

It is also possible to predict a dark matter halo's den- 
sity distribution if one assumes a mapping between a 
mass shell's initial radius, when the structure is linear, 
to its final average radius, which is some fraction of its 
turnaround radius. From this scheme, one can also infer 
a velocity dispersion, using the virial theorem. Unfortu- 
nately, this scheme does not give any information about 
the halo's velocity anisotropy. Our self-similar prescrip- 
tion discussed above, on the other hand, contains all of 
the velocity information. Hence, one can reconstruct the 
velocity anisotropy profile given the trajectory of a mass 
shell. The velocity anisotropy is significant since it de- 
scribes to what degree orbits are radial. Moreover, it can 
break degeneracies between n and w in our halo model. 
We will discuss the velocity structure of our halo model, 
including the pseudo-phase-space density profile, in more 
detail in Paper 2 of this series. There we will once again 
compare our halo predictions to the recent Aquarius sim- 
ulation results [6]. 

While the above self-similar prescription has its clear 
advantages, it's also unphysical since mass shells at 
turnaround are radially cold. The same tidal torque 
mechanisms which cause a tangential velocity dispersion 
[23 | . should also give rise to a radial velocity dispersion. 
For a more physical model, one would need to impose 
self-similarity to a phase space description of the halo 
and include sources of torque as diffusion terms in the 
Boltzmann equation. This will be the subject of Paper 3 
of this series. 

As we've shown, the way in which particles are torqued 
after turnaround (zu) influences the interior power law of 
the density profile. One way to source this change in 
angular momentum is through substructure that is as- 
pherically distributed throughout the halo. It is reason- 
able to assume that substructure dominated by baryons 
torque halo particles more strongly than substructure 
dominated by dark matter since baryons can achieve 
higher densities and hence are not tidally disrupted as 
easily. If this is the case, then torques sourced by baryons 
would result in a larger value of w than torques sourced 
by dark matter. According to the predictions of this sec- 
ondary infall model, this would lead to less cuspy profiles 
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(See Appendix [Cl for a more detailed discussion). There- 
fore a more thorough understanding of w coupled with 
this simplified model of halo formation could potentially 
shed light on the Cusp Core problem and thereby possi- 
bly bridge the gap between simulations and observations. 



Appendix B: Tidal Torque Theory 

We first derive eq. ([25)1 using cosmological linear per- 
turbation theory. Starting with the Zel'dovich approxi- 
mation [30l ] . we have: 
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where r is the physical radius, q is a Lagrangian coordi- 
nate, and 4>i which is time independent, is related to the 
Newtonian potential <£> through the following: 



Appendix A: Calculating n e ff 

The effective primordial power spectral index, n e g, re- 
lates our model parameter n to the halo mass M. The 
effective index is defined by: 



4:7rGp m Da 2 



(B2) 



where p m is the dark matter background density. The 
velocity, dr/dt, then is given by: 



_ dlna R 



where 



(Al) 



v(q,t) = H(t)r(q,t) - a(t)D(t)V4>(q) 



(B3) 



where H = dlna/dt and dots denote derivatives with 
respect to time. Therefore, to first order in 0: 



d 3 k 



P(k)W 2 R (k) 



(A2) 



and 



WR[k) = 



(kR) cos(fci?) 



(kR) 3 



(A3) 



The scale R is set by the top hat mass of the halo (M = 
47rp m o-R 3 /3), where p m o is the dark matter background 
density today. The power spectrum today, P{k), is given 
by m 



P(k) 



2k 2 



5Hq Q m 



P n (k)T 2 (k)D 2 (a = 1) (A4) 



where 



k 3 Pn(k) 
2tt 2 



= A^(fco) 



(A5) 



D is the linear growth factor normalized so that D/a — > 
1 as a — > 0, a is the scale factor, T(k) is the transfer 
function and we choose cosmological parameters derived 
from WMAP7: A^(fc ) = 2.441 x 1CT 9 , h = 0.704, fl m = 
0.272, n s = 0.963, with k = 0.002 Mpc -1 g|. We 
calculate T(k) using CMBFAST 0. 



(r - r ) X (v - v ) = -a 2 D(q - q ) x (V0 - V^o) 

(B4) 

where V0o = V0(qo)- Plugging this expression into eq. 
(|27| . taking an expectation value, defining x = q — q , 
rewriting in terms of the velocity perturbation variable 
\I/ = —D'Vfj) an d using index notation, we find: 



n 4 D 2 

<a 2 >= — — ' " 3 ~ 



D 2 



^-ijk^ilm I d Xp m CL XjXi X 



Vt 



< 



(* fc (0) - V k (x)) (* m (0) - * m (a;)) > c 

(B5) 



where e«/- is the Levi-Civita tensor and Vl is the La- 
grangian volume of the halo. We assume Vl is spherical 
with radius x max . Note that to linear order in 0, the 
Lagrangian coordinate x is equivalent to a comoving co- 
ordinate. The subscript C denotes an expectation value 
taken over constrained gaussian fields since we want to 
average only over peaks in the density field. 

Using the formalism developed in [47]], we choose, for 
simplicity, to use the zeroth and first order derivatives of 
the smoothed density field to constrain the velocity per- 
turbation. Our treatment and conventions are identical 
to that used in Appendix A of Ma and Bertschinger [48[ . 
For more details, please refer to this reference. We find: 
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x 2 ff{x) (£(x)-2fj(x)Y 



I X-iXj 



7(x) ??(0)?7(x) 



<* i (x)* i (0)> o =(^ 
fdr/{x) 7(3;) 0:77(0) d7j(^) 



v da; 



crj da; 



(*<(o)*i(o)) = (4-^r)« 



(B6) 

(B7) 
(B8) 



where: 



2 



I /" d 3 fc 



3 7 (2tt) 3 
d 3 k 



k- 2 P{k) 



0, 



' ~ l (2tt) 3 
1 /• d 3 fc 



77(2) 



1_ 3i (27T) 3 

d 3 fc 



k 2 P(k)W 2 R (k) 



P(k)W R (k) 



kx 



(2tt) 3 ' 

/d 3 k 
P(k)W R (k)j (kx) 

/H 3 k 
— P(k)k^ 3l {kx). 



(B9) 
(BIO) 
(Bll) 
(B12) 
(B13) 
(B14) 



W R {k) and P(k) are denned in eq. (|A"3"]l and (fAH) re- 
spectively and the spherical Bessel functions are jo 0*0 = 
x^ 1 sin a; and ji(x) = a; -2 (sin a; — a; cos a;). Eg. (|B7[) cor- 
rects an error in equation (A15) of Ref. |48j . 

Notice from eqs. (|B6|) through (|B8|) that the expres- 
sions separate into terms proportional to Sij and terms 
proportional to XiXj. The terms proportional to XiXj 
vanish in eq. (IB5|) because of the antisymmetry of the 
Levi-Civita tensors. Eq. (|B5|) then reduces to: 



(a 2 ) 



M 



t a 4 D 2 
1 D 2 



d*xp m a 3 x 2 f(x,R) (B15) 



where 



2 27(1) r7 2 (a;) 277(0)77(2;) ?7 2 (0) 



/Or, i?) = 24 



Last, defining u = x/x max , we find: 



(B16) 



Mr 

jyi max 



i f(ux max , R)di 



6a D Mx^^A (R) 



In the above A has units of Mpc and x max = i? since the 
scale of the galaxy today is equivalent to its lagrangian 
size to linear order in perturbation theory. Note that 
/ / D 2 is time independent. 

Now, we derive eq. (|5D|) . Consistency with the sec- 
ondary infall model demands that we assume Q m = 1. 
Equating the time dependence of eq. ([28]) and eq. ([29]) 
for an Einstein de-Sitter universe (D — a) at early times 
(r max oc t 2 / 3 ), we find p — In where p = 27 + 4. Given 
this relationship, we now equate eq. (|2"8")l to eq. (|2"91 and 
solve for £?. We find: 



B 



\y/2{l -2n)a 2 



'"la 



(B18) 



where we've used eqs. (flU)) and ([15} evaluated at early 
times to cancel the mass as well as the relationship 
r nlM = aa; max . Now we evaluate eq. (IB18I) at early times 
since tidal torque theory only applies when the halo is 
linear. To relate r max at some initial time (tj) to r ta 
today (to), we use the conservation of mass. 



4tt 



-PB(^>max(*i) 



47T 



-X(l)p s (t )r t 3 a (t ) (B19) 



Evaluating eq. (|B18j) at f j with the use of eq. (|B19j) and 
noting that r ta oc t@ , we find: 



B = ^ v /2(7-2n)X(l)(- 1 )/ 3 ^- 
3 rta (*o) 



(B20) 



As expected, the time dependence of B vanishes. Last, 
assuming r ta (io) = R, we reproduce eq. (|30[) . The quanti- 
ties n, R, A(R) are calculated in an VL m — 1 universe with 
the same background matter density and power spectrum 
of ACDM today. This ensures that the statistics, mass, 
and size of halos in both universes are equivalent today. 

As mentioned previously, eq. (f3T)|) overestimates the an- 
gular momentum of particles at turnaround by a factor of 
1.5 to 2.3. One potential source of error is to assume the 
lagrangian volume is spherical. Assuming an ellipsoidal 
lagrangian volume with axis ratios 1 : a : b, we find that 
B is at most reduced by 8% when 0.5 < a,b < 1. The 
discrepancy may be caused by not including higher or- 
der constraints on the smoothed density field. However, 
comparing B when calculated using zeroth and first order 
derivative constraints with B when calculated without 
using constraints results in only percent level differences; 
so it seems unlikely that constraints would have a signif- 
icant effect. 

N-body simulations use friends-of-friends group finders 
in order to identify halos and subhalos [49|, [50(. This 
algorithm, however, removes particles that are grouped 
to neighboring halos and hence neglects a contribution 
to er 2 . Trying to mimick this selection effect, we replaced 



(B17) eq. (jBTBl with 
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<- 2 > 



M 



1 D 2 



d 3 xe- x2/2R2 p m a 3 x 2 f(x,R) 



(B21) 

where R = (2/9 7 r) 1 / 6 i? ensures that the mass enclosed 
within the lagrangian volume is equal to the mass of the 
halo calculated by the simulation. However, calculating 
B in this manner leads to overestimating the angular 
momentum at turnaround by ~ 4, as opposed to ~ 2 
beforehand. Since f(x,R) is an increasing function of x 
(/ ~ x 1 - 34 near R for 10 12 M© halos), most of the con- 
tribution to A(R) comes from close to R. Therefore, 
while the gaussian cutoff decreases the contribution to 
B around R, it includes contributions beyond R, leading 
to a worse estimate. This highlights that B significantly 
depends on the outer parts of the halo. Hence, overesti- 
mating B by ~ 2 is reasonable. 

Last, note that the parameter B is set during the lin- 
ear regime. Assuming that the shell is dominated by 
substructure at turnaround, nonlinear interactions like 
dynamical friction and tidal stripping play an important 
role from the time of turnaround to the first pericenter 
passage [111, HH . As time goes on, these effects become 
less important since substructure in the shell becomes 
subdominant. Including these extra interactions should 
lead to smaller estimates of B at first pericenter passage 
and hence potentially explain our overestimate, but is 
beyond the scope of this work. 



Appendix C: Evolution After Turnaround 

In this Appendix, we use dimensional analysis in order 
to gain intuition about zu, & parameter that describes 
tidal torque after turnaround. First, consider the time 
derivative of L 2 , where L is the angular momentum per 
unit mass of a particle with radius r at time t. 



dL 2 
~dT 



2[r 2 (v ■ a) — (r ■ a)(r ■ v] 



(CI) 



In the above, a (v) is the acceleration (velocity) of the 
particle. We now decompose the acceleration vector into 
a radial (r) and tangential (t) component and use this 
basis to rewrite the velocity vector. 



a = a r r + a±t 

v — v r r + v t t + v p p 



(C2) 
(C3) 



The direction p is orthogonal to both f and t. Note 
that all basis vectors depend on position. Plugging in 
the above decomposed vectors into eq. (IC1|) . we find: 



dL 2 
~dT 



As expected, changes in L 2 are sourced by deviations 
from spherical symmetry that create nonzero a t . Now, 
imagine a spherically symmetric halo, roughly described 
by our self-similar infall model with vj > 0, with a clump 
of mass m in the shell at radius r-i. We assume that m 
is small enough so that it does not influence the radial 
equation of motion of the shell at r. We focus on in > 
since this is required in order for the density profile of 
a 10 12 M Q halo to be consistent with the NFW profile 
(Section MM- 

Next, consider averaging dL 2 / dt over an orbital period 
and over a spherical shell of radius r, in order to compare 
the change in angular momentum sourced by the clump 
at r 2 to our model's prescription for angular momentum 
evolution. For zu > 0, orbits are roughly circular at late 
times. Hence, we assume averaging over an orbital pe- 
riod is equivalent to evaluating the right hand side of eq. 
(|C4jl at roughly the apocenter radius of the shell. As 
described in Section |TT1 the orbital planes of all parti- 
cles in a shell at r are randomly aligned. Therefore we 
expect vt averaged over a sphere to vanish. However, if 
there exists an excess mass m, then all the particles will 
be pulled slightly in that direction, leading to a nonzero 
average. We therefore assume vt oc Pat where P, the 
orbital period of the particle, is taken to be a dynamical 
time (P oc p(r) -1 / 2 ). Last assuming r <C T2, we find: 



dL 2 
~df 



oc 



r 2 V /2 ' 



(C5) 



In order for the secondary infall halo model to be con- 
sistent, the right hand and left hand side must have the 
same time scaling. Assuming m oc i M , associating r and 
ri with their respective apocenters, and noting from eq. 
(1571) that r oc i 2ro , r 2 oc i 2ro , and p oc i" 6ro , we find: 



^=3(1 + 2^) 



(C6) 



2r v t a t . 



(C4) 



The same scaling relationship holds for r 3> r-z. Hence 
one could imagine substructure in the shell at r sourc- 
ing a change in L 2 of the shell at r 2 and substructure 
in the shell at r 2 sourcing a change in L 2 of the shell at 
r. Therefore, a hierarchy of substructure non-spherically 
distributed, which is subdominant to the monopole con- 
tribution of the halo, would result in a halo roughly con- 
sistent with our described secondary infall model. 

Eq. (|C6I) . which is only valid for p > —1/2 since we as- 
sumed w > 0, together with eq. (157)) relates the steepness 
of the inner density profile to the mass loss rate of sub- 
structure. If the clump does not lose mass (p = 0), then 
zu = 1/3 implies p oc r . If the clump loses mass {p < 0), 
eq. (|5Tj) predicts steeper density profiles. Substructure 
dominated by baryons will lose less mass than substruc- 
ture dominated by dark matter, since baryons clump 
more easily and hence have higher densities. Therefore, 
according to the above analysis, pure dark matter simu- 
lations should have steeper density profiles than galaxies 
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which include baryons. This is expected since baryons 
stir particles around more efficiently, causing larger peri- 



centers and less dense interiors. A more thorough treat- 
ment that involves constraining w with simulations is 
beyond the scope of this paper. 



[1] 

[2] 
[3] 
[4] 
[5] 

[6] 



[7] 
[8] 
[9] 

[10] 

[11] 
[12] 

[13] 
[14] 
[15] 
[16] 



J. E. Gunn and J. R. Gott, III, Astrophys. J. 176, 1 
(Aug. 1972). 

J. R. Gott, III,|Astrophys J.| 201, 296 (Oct. 1975). 
J. E. Gunn, [Astrophys. J.|218 , 592 (Dec. 1977). 

E. Bertschinger, |Astrophys. J.| 58, 39 (May 1985). 

J. A. Fillmore and P. Goldreich, | Astrophys. J.| | 281, 1 
(Jun. 1984). 

J. F. Navarro, A. Ludlow, V. Springel, J. Wang, M. Vo- 
gelsberger, S. D. M. White, A. Jenkins, C. S. Frenk, 
and A. Helmi. IMon. Not. R Astron. Socl 402. 21 (Feb. 
2010), larXiv:0810.1522l 

A. W. Graham, D. Merritt, B. Moore, J. Diemand, 
and B. Terzic, lAlitron. J.I 132, 2701 (Dec. 2006), 
|arXiv:astro-ph/ 0608613 

J. Diemand, B. Moore, and J. Stadel, 
IMon. Not. R. Astron. Socl 353, 624 (Sep. 2004), 
|arXiv:astro-ph /0402267[ 

W. J. G. de Blok, in Revista Mexicana de Astronomia 
y Astrofisica Conference Series, Revista Mexicana de 
Astronomia y Astrofisica, vol. 27, Vol. 17, edited by 
V. Avila-Reese, C. Firmani, C. S. Frenk, & C. Allen 
(2003) pp. 17-18. 

G. Gentile, P. Salucci, U. Klein, D. Vergani, and 
P. Ka lberla. IMon. Not. R. Astron. Sod 351, 903 (Jul. 
2004), |arXiv:astr o-ph/040 3T54| 

P. Salucci, A. Lapi, C. Toiiini, G. Gentile, I. Yegorova, 
and U. Klein, IMon. Not. R. Astron. Sod 378, 41 (Jun. 
2007), |arXiv:astro-ph/0703115j 

F. Donato, G. Gentile, P. Salucci, C. Frigerio Martins, 
M. I. Wilkinson, G. Gilmore, E. K. Grebel, A. Koch, and 
R. Wy se, IMon. Not. R. Astron. Soc.l 397, 1169 (Aug. 
2009), |arXiv:0904.4054 [astro-ph.CO] | 

B. S. Ryden and J. E. Gunn, |Astrophys. J.| 318, 15 (Jul. 
1987). ' 

A. Nusser, IMbn. Not. R. Astron. Soc.l 325, 1397 (Aug. 



2001), |arXiv:astr o-ph/0008217 
N. Hiotelis, |Astron. Astrophys.| 382, 84 (Jan. 2002), 
|arXiv:astro-ph/0111324"r 

L. L. R. Williams, A. Babul, and J. J. Dal- 
canton, |Astrophys."jT| 604, 18 (Mar. 2004), 
|arXiv:astro-ph/0312002[ 



[17] 
[18] 
[19] 

[20] M. Le Delliou and R. N. H enriksen , |Astron. Astrophys.| 
[21] 



P. Sikivie, I. I. Tkach ev, and Y. Wang,|Phys. Rev. D| 56, 

1863 (Aug. 1997 ) ,|arXiv:astro-ph/9609022| 

A. Del Popolo, |Astrophys. J.| 698, 2093 (Jun. 2009), 

larXiv:0906.4447l 

S. D. M. White and D. Zaritsky, |Astrophys.~T| 394, 1 
(Jul. 1992). 



408, 27 (Sep. 2003), arXiv:astro-ph/0307046 
Y. Ascasibar, G. Yepes, S. Gottlober, and V. Miiller, 
IMon. Not. R. Astron. goal 352, 1109 (Aug 



2004), 

|arXiv:astro-ph/031222T] 

[22] Y. Hoffman and J. Shaham, [Astrophys. J.| 297, 16 (Oct. 
1985). 

[23] F. Hoyle, in Problems of Cosmical Aerodynamics (1951) 
pp. 195-197. 



[24 

[25 
[26 

[27 
[28 
[29 
[30 
[31 
[32 

[33 

[34 

[35 

[36 

[37 

[38; 

[39 

[40 

[41 
[42 
[43 
[44 
[45 



V. Springel, 
(May 2007), 



E. Hayashi, J. F. Navarro, and 
IMon. Not. R. Astron."5ocT1 377, 50 

|arXiv: astro-ph /06 12327] 

J. Barnes and G. Efstathiou, Astrop hys. J.| 319, 575 
(Aug. 1987). 

M. Boylan-Kolchin, V. Springel, S. D. M. White, and 
A. Je nkins, IMon. Not. R. Astron. Soc.l 406, 896 (Jun. 
2010) , larXiv:091 1.44841 

P. J. E. Peebles, [Astrophys. J.I 155, 393 (Feb. 1969). 
S. D. M. White, [Xstrophys. J.[ 286, 38 (Nov. 1984). 
A. G. Doroshkevich, Astrofizika 6, 581 (1970). 
Ya. B. Zel'dovich, Astron. Astrop hys. 5, 84 (Mar. 1970). 
S. Chandrasekhar, |Astrophys. J.| 97, 255 (Mar. 1943). 
O. E. Gerhard and J. Binney, Mon. Not. R. Astron. Soc. 
216, 467 (Sep. 1985). 

D. Merritt and G. D . Quinla n, |Astrophys. J7| 498, 625 
(May 1998), |arXiv:astro-ph/9709106[ 

F. Cruz, H. Velazquez, and H. Aceves, Revista Mexicana 
de Astronomia y Astrofisica 43, 95 (Apr. 2007). 

A. J. Kalnajs, in Dynamics of Disc Galaxies, edited by 

B. Sundelius (1 991) pp. 323 -+. 

W. Dehnen, lAstronTjl 119, 800 (Feb. 2000), 
|arXiv:astro^ph/991116l1 

M. Milosavlj evic and D. Merritt, |Astrophys. J.| 596, 860 
(Oct. 2003), [arXiv:astro-ph/0212459| 



A. Sesana, F. Haardt, and P. Madau,|AstrophysT j, 660, 
546 (May 2007), [arXiv:astro-ph/0612265[ 

J. F. Navarro, C. S. Frenk; and S. D. M. 
White, | Astrophys. J.I 462, 563 (May 1996), 
|arXiv:astro-ph/9508025] 

B. Moore, T. Quinn, F. Governato, J. Stadel, and 
G. La ke, iMoii. Not. R. Astron. Soc.l 310, 1147 (Dec. 

1999), |arXiv:astr o-ph/990 3164[ 

M. Ricotti, A. Pontze n, and M. Viel, |Astrophys. J. Lett.| 
663, L53 (Jul. 2007h larXlv:0706.0856l " 

R. Cen, F. Dong, P. Bode, and J. P. Ostriker, ArXiv As- 
trophysics e-print s(Mar. 2004), |arXiv:astro-ph/0403352"| 
A. Del Popolo, IMon. Not. R. Astron. Soc.l 1312T5ep. 
2010). 

M. Takada, E. Komatsu, and T. Futamase, Phys. Rev. D 
73, 083520 (Apr. 2006), [arXiv:astro-ph/0512374^ 



E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett, 
B. Gold, G. Hinshaw, N. Jarosik, D. Larson, M. R. 
Nolta, L. Page, D. N. Spergel, M. Halpern, R. S. Hill, 
A. Kogut, M. Limon, S. S. Meyer, N. Odegard, G. S. 
Tucker, J. L. Weiland, E. Wollac k, and E . L. Wright, 
ArXiv e-prints(Jan. 2010), i arXiv: 1001. 45381 
[46] U. Seljak a nd M. Zaldarriaga, | Astrophys. T\ 469, 437 

(Oct. 1996), [arXiy:astro-ph/9603033| 
[47] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, 

|Astrophys."JJ i 304, 15 (May 1986). 

[48] C. Ma and E. B ertsching er, |Astrophys. J.| 612, 28 (Sep. 

2004) , |arXiv:astro-ph/0311049T ~ 
[49] M. Davis, G. E fstathiou, C. S. Frenk, and S. D. M. White, 

|Astrophys. J . 292, 371 (May 1985). 
[50] V. Springel, J. Wang, M. Vogelsberger, A. Ludlow, 



18 



A. Jenkins, A. Helmi, J. F. Navarro, C. S. Frenk, and 
S. D. M. White, IMon. Not. R Astron. Soc.l 391, 1685 
(Dec. 2008), arXiv:0809.0898] 
[51] J. Gan, X. Kang, F. C. van den Bosch, and J. Hou, ArXiv 



e-prints(Jun. 2010). larXTv:1007.0023l 
[52] A. R. Zentner, A. A. Berlind, J. S. Bullock, A. V. 
Kravtsov, and R. H. Wechsler, Astrophys. J. 624, 505 
(May 2005), |arXiv:astro-ph/0411586| ~ 



